subroutine ReadMol_init()
use molecule
implicit none
   !-------------------------------------------------------------------------
   ! File format
   ! line 1: molecular name
   ! line 2: # of atoms in mol, max # of neighbors for each atom
   ! line 3: commen line
   ! line 4 -> natom+3: id EName x y z q neigh_id1 neigh_id2 ...
   !-------------------------------------------------------------------------
   character (len=256) :: molfile
   character (len=512) :: oneline, append
   integer             :: ioerr, id, i, idum
   logical             :: fexst
   !-------------------------------------------------------------------------
   write(*, '(/,10x,"Please input the filename that contains the mol info: ", $)')
   read(*, *, iostat=ioerr) molfile
   if (ioerr.ne.0) stop
   !
   inquire(file=molfile, exist=fexst)
   if (.not.fexst) then
      write(*,'(10x,"File:", A, " not found!")') molfile
      stop
   endif
   open(10, file=molfile, action='read')
   read(10, *, iostat=ioerr) molname
   if (ioerr.ne.0) stop
   read(10, *, iostat=ioerr) natom, neighmax
   if (ioerr.ne.0.or.natom.lt.1.or.neighmax.lt.1) then
      write(*,'(10x,"Error while reading natom and neighmax from ",A)') molfile
      stop
   endif
   read(10,*)
   !
   if ( allocated(neigh)  ) deallocate( neigh   )
   if ( allocated(element)) deallocate( element )
   if ( allocated(molpos) ) deallocate( molpos  )
   if ( allocated(chg)    ) deallocate( chg     )
   allocate( element(natom), neigh(neighmax, natom), molpos(3,natom), chg(natom))
   !
   append = ""
   do i = 1, neighmax
      append = trim(append)//" 0"
   enddo
   do i = 1, natom
      read(10, '(A)', iostat=ioerr) oneline
      if (ioerr.ne.0) stop
      oneline = trim(oneline)//trim(append)
      read(oneline, *, iostat=ioerr) id
      read(oneline, *, iostat=ioerr) idum, element(id), molpos(:,id), chg(id), neigh(:,id)
      if (ioerr.ne.0) stop
   enddo
   close(10)
   !-------------------------------------------------------------------------
return
end subroutine ReadMol_init
